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ABSTRACT 

The process referred to as "semi-convection" in astrophysics and "double- 
diffusive convection in the diffusive regime'' in Earth and planetary sciences, 
occurs in stellar and planetary interiors in regions which are stable according to 
the Ledoux criterion but unstable according to the Schwarzschild criterion. In 
this series of papers, we analyze the results of an extensive suite of 3D numerical 
simulations of the process, and ultimately propose a new ID prescription for heat 
and compositional transport in this regime which can be used in stellar or plane- 
tary structure and evolution models. In a preliminary study of the phenomenon, 
Rosenblum et al. (2011) showed that, after saturation of the primary instability, 
a system can evolve in one of two possible ways: the induced turbulence either 
remains homogeneous, with very weak transport properties, or transitions into 
a thermo-compositional staircase where the transport rate is much larger (albeit 
still smaller than in standard convection). 

In this paper, we show that this dichotomous behavior is a robust property 
of semi-convection across a wide region of parameter space. We propose a simple 
semi-analytical criterion to determine whether layer formation is expected or 
not, and at what rate it proceeds, as a function of the background stratification 
and of the diffusion parameters (viscosity, thermal diffusivity and compositional 
diffusivity) only. The theoretical criterion matches the outcome of our numerical 
simulations very adequately in the numerically accessible "planetary" parameter 
regime, and can easily be extrapolated to the stellar parameter regime. 

Subsequent papers will address more specifically the question of quantifying 
transport in the layered case and in the non-layered case. 
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Introduction 



1.1. The physics of semi-convection: double-diffusive convection 



The concept of "semi-convection", first introduced by Schwarzschild & Harm (1958), is 



often invoked in a number of rather different situations ( Merryfield||1995 ) which nevertheless 
have one point in common: they occur in regions which are stable to the Ledoux criterion, 
but unstable to the Schwarzschild criterion. Mathematically speaking, this condition can be 
expressed as 

'd\nT\ fd\nT\ fdln^i 



< 



< 



equivalently < V — Vad < V, 



(1) 



dlnp J \d\np J \d\np 

where T, p and fi are the temperature, gas pressure, and mean molecular weight, and where 
the subscript "ad" denotes a derivative at constant entropy. Physically speaking, ([T]) de- 
scribes regions which are thermally unstably stratified, but where standard convection is 
suppressed by the presence of significant compositional gradients. 

The first linear analysis of the stability of "semi-convective" regions was presented by 



Walin (1964) in the oceanographic context, and later by Kato (1966) in the astrophysical 



context. They both showed that a semi-convective region is hydrodynamically unstable, 
but to a much more gentle instability than the one associated with standard convection 
because it relies on doubly-diffusive processes to grow. It is in fact one of two forms of so- 
called "double-diffusive convection" (the other being fingering convection, otherwise known 
as "thermohaline convection" in astrophysics), and is often referred to as "double-diffusive 
convection in the diffusive regime" in physical oceanography. For the sake of clarity and 
brevity, we will refer to the phenomenon as "diffusive convection" in this series of papers. 



As reviewed by Rosenblum et al. (2011), diffusive convection is principally controlled 



by three non-dimensional parameters. The first two characterize the nature of the fluid 
considered, and are the Prandtl number Pr and the diffusivity ratio r, 

Pr = — , r = ^ , 

where u, kj- and are the viscosity, thermal diffusivity and compositional diffusivity re- 
spectively. For reference, note that Pr and r are very roughly of the order of 10~^ in giant 
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planet interiors and 10 ^ in stellar interiors, and that r is usually somewhat smaller than 
Pr. 



The third parameter is the inverse density ratio, defined as 



V- V 



ad 



(3) 



which measures the relative importance of the destabilizing thermal stratification compared 



with the stabilizing compositional one. A semi-convective region is unstable (Walin 1964 

Pr 



Kato 1966) if 



1< i?n ^ < 



1 



Pr 



(4) 



whereas regions with Rq^ < 1 are unstable to overturning convection and those with R^^ > 
R~^ are absolutely stable. Within this parameter range, it can be shown that the growth rate 
of the linear modes is complex. Furthermore, in the low Prandtl number regime characteristic 
of stellar and planetary interiors, the real part of the growth rate is proportional to the square 
root of the Prandtl number times the Briint-Vaisala frequency (see Appendix A). In the same 
limit, the typical lengthscale of the unstable modes is a thermal diffusion scale. 

Linear stability is unfortunately of rather limited utility, in particular when it comes to 
estimating the mixing rates induced by the diffusive convection. Experiments - laboratory or 
numerical - are the only way forward. Since no terrestrial fluid exists with similar values of 
Pr and r, it is tempting to use experimental measurements of heat and compositional fluxes 
in different parameter regimes, in particular laboratory experiments in the heat-salt system 



relevant for physical oceanography (Linden & Shirtcliffe 1978), and extrapolate them to the 



astrophysical case ( |Stevenson|p82| [Guillot et al.|[2004| . However one must be very cautious 
in doing so since the Prandtl number of water is typically 4-7, depending on temperature, 
whereas the Prandtl number in stellar and planetary interiors is much lower than one. Since 
the typical scale of the instability is a thermal diffusion scale, double-diffusive mixing is a 
mostly-laminar process at high Pr and a turbulent one at low Pr. There is no reason to 
expect that the laminar scalings should apply to the turbulent case. 

Nevertheless, from a qualitative point of view, one of the most interesting results from 
laboratory ( Turner fc Stommel||1964 ) and fleld experiments ( Timmermans et al.||2008 ) in salt 
water is the fact that diffusive convection has a tendency to form thermohaline staircases, 
i.e. well-deflned mixed layers separated by thin, very strongly stratifled, and essentially 
diffusive interfaces. It is often thought that the necessarily weak transport through these 
interfaces is what controls and limits the efficiency of transport by "layered" convection. Such 



considerations have led Spruit (1992) and Chabrier & Baraffe (2007) for example to propose 



theories for heat and compositional transport in astrophysics, which rely on assumptions 
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about the layer heights and the interface thicknesses. However, it is important to remember 
that, until very recently, layered convection had never been demonstrated to exist at low 
Prandtl number. As a result, these theories have, by and large, remained un-tested (see 



Rosenblum et al. (2011) for a review of prior numerical work). 



1.2. Recent numerical and theoretical results 

Recent 3D numerical simulations have finally shed some light on the subject of transport 



by diffusive convection. Rosenblum et al. (2011) ran an exhaustive numerical study of the 



phenomenon for fixed Prandtl number and diffusivity ratio Pr= r = 0.3, and with the 
inverse density ratio Rq^ spanning the entire instability range (j4]). While still far from 
any astrophysically-relevant regime, the selected values of Pr and r in these simulations 



were below unity and therefore in the right "region" of parameter space. Rosenblum et al. 



(2011) found that the instability grows as expected and that the early behavior of the fluid 
can be satisfactorily explained by considering the fastest-growing modes of instability only. 
However, they discovered that two very different regimes of diffusive convection are possible 
after saturation of what we will refer to as the "primary" instability, depending on the value 
of the inverse density ratio. The various regimes are illustrated in Figure [1} 



For large enough Rq , i.e. for more "stable" stratifications, Rosenblum et al. (2011) 
showed that the system settles into a homogeneous, statistically stationary, weakly turbu- 
lent state. The turbulence is dominated by internal gravity waveband mixing occurs princi- 
pally via wave-breaking. Heat and compositional transport are fairly inefficient, and depend 
sensitively on the inverse density ratio. For low enough i?g ^ on the other hand, which cor- 
responds to systems closer to the Ledoux-stability criterion, they observed the spontaneous 
emergence of thermo-compositional staircases after a short adjustment period. These stair- 
cases take the form of vigorously convective layers, which are thermally and compositionally 
well-mixed, and separated by fairly sharp interfaces, as observed in the oceanographic case. 
The interfaces, however, are far from merely diffusive and are instead very dynamic, and 
often pierced by strong localized updrafts and downdrafts. Later on, the layers are observed 
to merge, and each merger is accompanied by a significant increase in the overall transport 



across the staircase. Rosenblum et al. (2011) found that transport in that regime depends 



sensitively on the mean layer height, and proposed preliminary scalings to quantify it. The 
latter remain to be verified across a wider region of parameter space. 



^Recall that the background stratification is stably stratified in terms of the density and therefore supports 
standard gravity waves which oscillate with the buoyancy frequency. 
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Fig. 1. — Illustration of the various regimes of diffusive convection. In systems without com- 
positional gradients, the Schwarzschild criterion marks the stability boundary between overturn- 
ing convection and absolute stability. In the presence of a stable compositional gradient, diffu- 
sive convection occurs for Rq^ between 1 (which corresponds to the Ledoux-stability limit) and 
= (Pr + l)/(Pr-Fr). Within this range, two possibilities arise: for Rq^ €^ [1,R^^], spontaneous 
transition into layered convection is observed, while for Rq^ G I^l^' -^cT^Ji ^he system remains in a 
state of weak oscillatory convection. Note that both R~^^ and R^^ depend on Pr and r. 
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In fact, insight into the reason for this dichotomous "layered vs. non- layered" con- 
vection, can be gained from recent oceanographic studies of fingering convection, a related 
double- diffusive instability of thermally stable fluids which are destabilized by adverse com- 



positional gradients (Stern 1960). Such conditions are found in the tropical thermocline for 



example, where surface heating and evaporation continually warm up the upper layers of wa- 
ter, and increase its salt concentration. Crucially, thermohaline staircases are also commonly 
found in fingering regions of the ocean dSchmitt et aL]|2005D . |Radko| ( |2003j ) studied their for- 
mation in this regime, and showed that they can naturally emerge as a secondary large-scale 
"mean-field" instability of the system. More precisely, he showed that horizontally invariant 
but vertically sinusoidal density perturbations grow exponentially out of the homogeneous, 
small-scale fingering convection, and eventually overturn into a regularly-spaced staircase. 



His theory was later validated by Stellmach et al. (2011) via three-dimensional numerical 
simulations. 

A crucial result of Radko's theory is his identification of a necessary and sufficient 
condition for the layer-forming instability, namely that the turbulent flux ratio 7turb5 defined 
as the ratio of the turbulent buoyancy flux due to heat transport, to the turbulent buoyancy 
flux due to salt transport, should be a decreasing function of the density ratio Rq. He 
thus named the instability "the 7— instability" . The fact that 7turb in salt water decreases 
with Rq for low density ratios, then increases again for higher density ratios, explains why 
thermohaline staircases in the tropical ocean are only found in regions with low enough Rq. 



Recently, Rosenblum et al. (2011) showed that Radko's 7— instability theory can very 



easily and naturally be extended to explain the emergence of staircases in their own simu- 
lations of diffusive convection. The equivalent condition for instability is that the total flux 
ratio 7tot (i-e. the ratio of the total buoyancy fluxes, diffusive plus advective, of heat to 
composition respectively), should be a decreasing function of Rq. Since the inverse density 
ratio Rq^ is a more convenient parameter in diffusive convection, an equivalent sufficient 
condition for instability is that 7^"?; should be a decreasing function of Rq^- For complete- 
ness, the 7— instability theory is rederived and discussed in Section |3] Rosenblum et al. 
(2011) found through their systematic exploration of the instability range that 7^"^ has a 

0.3. This explains why layers are seen to 



minimum at about Rj^^ 



1.4 when Pr = r 



form for Rq^ < 1.4 in their simulations (at these values of Pr and r) but not for larger 



Rq . Furthermore, in the simulations which do lead to layering, Rosenblum et al. (2011) 



found that theory and numerical experiments agree remarkably well on the growth rate of 
the 7— instability. Their preliminary study thus suggested that, in order to know under 
which conditions layer formation is possible in stars and giant planets, one simply needs to 
determine if and when 7^"^; decreases with Rq^, for a given parameter pair (Pr, r). 
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1.3. Work outline 



The findings of Rosenblum et al. (2011 ) lay a very clear path towards creating a practical 



model for transport by diffusive convection (semi-convection) in astrophysics: 

1. Model the function ■j^oti^o^j ''") (from numerical simulations and/or theoretical cal- 
culations) to determine if and when layered convection is expected, and verify whether 
the 7-instability predictions continue to hold at lower Pr and r. 

2. Characterize transport by layered convection, and in particular, its dependence on 
layer height, Pr, r and Rq^- 

3. Characterize transport by homogeneous diffusive convection (i.e. in the absence of 
layers) . 

The first step is addressed in this paper, while steps 2 and 3 are deferred to subsequent 
publications in the series. 

In the present paper, we therefore focus our efforts on a precise determination of the 
region of parameter space where layered convection is expected to occur. We do so using 
a combination of numerical simulations and theory. We discuss the numerical model and 
present typical results in Section |2j We review the 7— instability theory in Section [3j In 
Section [i] we outline the methodology used for extracting the value of the flux ratio 7^"^ 
from simulations at numerically-accessible parameters, and present our results. In Section 
[5] we present a simple semi-analytical theory which enables us to estimate 7^" J for any set 
of parameters, and compare it with our numerical results. In Section [6} we show that the 
predicted growth rates from the 7— instability theory match the results of our numerical 
simulations very well. Finally, we conclude in Section [7) 



2. Mathematical model and typical solutions 
2.1. Mathematical model 



As discussed by Rosenblum et al. (2011), the typical lengthscale of the fastest unstable 
modes is of the order of meters to hundreds of meters at most in the parameter regimes typical 
of stellar and planetary interiors. They found that the first layers to form are quite thin, 
spanning no more than a few of these fastest-growing wavelengths. This justifies studying 
diffusive convection (at least, in the early stages of layer formation and evolution), as a local 
rather than a global process. 
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We consider a local Cartesian domain of size {L^, Ly, L^), where gravity defines the 
vertical direction: g = —ge^. The small domain size permits the use of the Boussinesq 
approximation ( Spiegel fc Veronis|1960 ) to the governing equations, which are then expressed 



as 



Vu 



dT 
~dt 



+ u ■ VT + (To, - T, 



ad- 
Oz , 



W 



+ U ■ VyU + llQzW 



, 

ktV^T 
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+ u ■ Vu 



Po 



(5) 



The first of these equations is the continuity equation, where u = {u,v,w) is the velocity 
field. The temperature and chemical composition fields (the latter is represented here for 
example as the mean molecular weight of the fluid) are expressed as the sum of a linear 
background profile {zTq^ and zfioz) plus triply-periodic perturbations T and fi. The thermal 
energy equation has been re-written as an advection-diffusion equation for T, with kt being 
the thermal diffusivity. The additional term —wTq^ is present in compressible fluids but not 
in incompressible ones, and represents the temperature change due to adiabatic expansion 



(Spiegel & Veronis 1960). Another advection-diffusion equation models the evolution of 



the mean-molecular weight perturbation /x, with the corresponding diffusivity. The last 
equation in ^ is the momentum equation; in the Boussinesq approximation, the density 
perturbation about hydrostatic equilibrium, p, appears in the buoyancy term only, and is 
linearly related to the temperature and mean molecular weight perturbations as 



Po 



-aT + 



(6) 



where po is the (constant) mean density of the region considered, and a and (3 are the 
coefficients of thermal expansion and compositional contraction respectively. The pressure 
perturbation is denoted as p, and u is the viscosity. All the perturbations satisfy triply- 
periodic boundary conditions. 



q{x, y, z, t) = q{x + L^., y, z, t) = q{x, y + Ly, z, t) = q{x, y,z + L^, t) , 
where q G {u, T, This setup minimizes the effects of boundaries on the system. 



(7) 



Using the following standard non-dimensionalization (jRosenblum et al. 2011|): 



[l]=d 



ag\T,,-T-f\ 



1/4 



- 9 - 



[t] = d^KT , 

[T] = t/|ro.-Tof| , 

M = («//3)|To,-Tof|rf, 

the governing equations can be re-written as 

u-Vu^ = -Vp + (f - /i)e^ + V^u 

u-Vf-w = V^T , 
+ u ■ V/i - Ro^w = rV^/i , 



1. 

df 
~dt 

di 







(9) 



where quantities with tildes are dimensionless. The three parameters discussed in Section 



naturally appear, namely the Prandtl number Pr, the diffusivity ratio r, as well as the 
inverse density ratio -Rq ^, see equations pi) and O). In the notations used here, we also have 



(10) 



p-l _ /^/^Oz 



It is interesting and important to note that this non-dimensional model now only knows 
about the superadiabatic temperature gradient Tq^ — T^^ rather than about Tq^ and Tgf 
individually. As such, any two real physical systems with the same superadiabaticity, the 
same density ratio, and the same values of Pr and r, will lead to the same non-dimensional set 
of equations even if their background temperature gradients are different. This degeneracy 
in the parameters will be discussed in more detail in Section |3| 



All simulations presented in this work were obtained using the PADDI code (see Traxler 



et al.| ( |2011b| ; |Rosenblum et al.| ( |201l| )), which solves in a cubic domain of size (lOOd)" 
subject to boundary conditions ([T]) using pseudo-spectral DNS. The selection of the domain 
size is discussed in Section HI 



2.2. Typical results 

Here we present the results of two selected simulations, which illustrate the behavior 



of diffusive convection in the planetary parameter regime. As mentioned in Section 1.2 



previous work at Pr= r = 0.3 showed that the evolution of the system after saturation of 
the primary double-diffusive instability can either result in layer formation or not, depending 
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on the value of the inverse density ratio Rq^- We confirm that this is still true at lower Pr 
and T. This is shown in Figure |2| which illustrates the two regimes for Pr = r = 0.03: the 
layered case, using R-^ = 1.5 (top row), and the non-layered case, using Rq^ = 5 (bottom 
row). In each case, we show on the left a snapshot of the simulation at a particular time t. On 
the right we show the temporal evolution of the non-dimensional thermal and compositional 
fluxes (wT) and (wp,), where (■) denotes a spatial average over the entire computational 
domain. 

In the case with low Rq^ (top row), the basic instability rapidly saturates (around 
t = 500 in non-dimensional time units), then transitions into a layered state, around t = 
1200. Three easily identifiable layers initially appear, which then merge into two (at about 
t = 1500) then one (at about t = 1800). The snapshot on the top left of Figure |2] was 
taken at t = 1760, and shows the non-dimensional perturbation in the concentration field in 
the 2-layered state. The fiuxes clearly increase in a stepwise manner, first when layers form 
and then at each merger. This kind of behavior was already illustrated and discussed by 



Rosenblum et al. (2011) (see their Figure 7). 

In the case with high Rq^, the basic instability also grows (although more slowly, as 
expected from linear stability) and eventually saturates around t = 3000. However, layers 



never form. Instead, what follows saturation is what Rosenblum et al. (2011) described as 
being a homogeneous, weakly convective phase with fairly inefficient transport properties. 
The snapshot on the bottom left of Figure [2] shows the non-dimensional perturbation in 
the concentration field at t = 4600. Note the small amplitude of the perturbations, by 
comparison with the total compositional contrast across the domain (A/z = 500 here). Upon 
closer inspection, we find that the small-scale oscillatory structures that are characteristic of 
the homogeneous phase intermittently give way to somewhat larger-scale and more coherent 
gravity waves (e.g. here for 4700 <t< 5200 and t > 5500; see also Figure |3|. When this is 
the case the amplitude of the wave-induced oscillation in the fiuxes dramatically increases, 
and the mean wave-induced transport also increases, although remains much lower than 
in the layered phase. The reason for the emergence of larger-scale waves and their self- 
organization remains to be determined, but the phenomenon is fairly ubiquitous at high R^^ 
(see below). This effect will be studied in a subsequent paper. 



Figure [3] shows the evolution of the turbulent heat flux for parameter pairs (Pr, r) with 
Pr = r, for selected Rq^ ranging from values close to overturning instability (left column), 
through intermediate values (middle column) to values close to marginal stability (right 
column). The plot clearly illustrates the following trends. Simulations with the lowest 
values of Rq^ lead to very rapid layer formation, while those with slightly larger values of 
Rq^ can stay in a state of homogeneous diffusive convection for a very long time before 
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Fig. 2. — Example of simulation results for Pr = r = 0.03, for Rq^ = 1.5 (top row) and Rq^ = 5 
(bottom row). The figures on the left are snapshots of the compositional perturbation field, at 



t = 1760 for the R, 







1.5 case and t 



4600 for ^ 



5 case. Note the vast difference in the 



amplitude of the perturbations for the two cases: for reference, the total compositional contrast 
across the domain is = 150 for Rq^ = 1.5 and A^u = 500 for Rq^ = 5. As a result, the 
density profile has local inversions in the low R^^ case (i.e. "layers"), but remains very close to the 
background state in the high Rq^ case. The figures on the right show the corresponding temporal 
evolution of the non-dimensional turbulent fluxes {wT) and {wji). Note the stepwise increase in 
the layered case, with layer formation and each subsequent merger. 
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layers emerge (see the case of Pr = r = 0.1, i?Q = 1.5 for example). At intermediate 
values of -Rq ^' layers never form. The system remains in a state of homogeneous diffusive 
convection, and occasionally exhibits intermittent gravity-wave-dominated phases similar to 
the one described earlier. Finally, in runs with larger values of R^^ closer to marginally 
stability we always see that the wave-dominated phase begins very quickly after saturation 
of the primary instability. These various types of behavior need to be kept in mind when 
analyzing the data to extract their mean transport properties, as described in Section |4} 



3. The 7— instability revisited for case Vad 7^ 

In this Section we rederive the 7— instability theory for the sake of completeness, and 
to correct a slight inconsistency in nomenclature discovered in the work of [Rosenblum 



et al. (2011). While their derivation is technically correctj^ it requires a slight physical 



re-interpretation of the quantities they define as "thermal Nusselt number" and "total buoy- 
ancy flux ratio" in order to be fully consistent when Vad 7^ 0, as shown below. 



The 7— instability theory, first proposed by Radko (2003) in the context of fingering 
convection in the ocean, is a mean-field theory that describes the development of secondary 
instabilities in fully-developed double-diffusive convection (the theory is in fact valid both 
in the diffusive regime and in the fingering regime). The theory assumes that the system is 
already in a homogeneous and quasi-steady turbulent state and studies its evolution when 
subject to perturbations on scales much larger than the turbulent eddies. Accordingly, we 
begin by averaging the governing non-dimensional equations (|9]) over all small lengthscales 
and fast timescales, and study the evolution of the large-scale, more slowly evolving mean 
fields. 

An emerging staircase is a horizontally invariant structure with no mean flow. If we 
ignore the momentum equation, and neglect mean flows as well as horizontal derivatives, the 
averaged non-dimensional thermal and compositional advection-diffusion equations become: 



dt dz 



tot 



dt dz 



where • denotes a spatio-temporal average over small-scales and short timescales. Note that 
the vertical fluxes and include a diffusive and a turbulent component. The goal is 



^See associated Erratum for correction of a typo. 
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Fig. 3. — Temporal evolution of the turbulent heat flux (wf) for parameter pairs (Pr,r) with 
Pr = r = 0.3, 0.1, 0.03 and 0.01 respectively from top to bottom in each column. In each plot, 
the X— axis represents the non-dimensional time t, and the values of Rq^ corresponding to each 
run represented are indicated. The left-column only shows runs which arc found to transition into 
layers. The middle and right columns show runs at intermediate and high values of Rq^ respectively. 
Runs at the highest values of Rq^ are often immediately dominated by large-scale gravity waves. 
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to express them in terms of large-scale fields only and thus close the system of equations, so 
that the latter can be solved for the evolution of T{z,t) and Jl{z,t). 

It is important to note for the upcoming discussion that there is, at this point, some 
degree of flexibility in the definition of these two fluxes: one can add or subtract any constant 
to and without changing the expression dF^^^/dz. In the original derivation of the 
7— instability, the fluxes are thus taken to be the total non-dimensional heat and composi- 
tional fluxes through the system, including the diffusion of the background fields Tq^ and 
fioz- When expressed non-dimensionally. 



Kt\I0z - -l-oz 



nAa/(3)\To.-T,i\ + ^^^^ ' 

where the subscript z denotes a derivative with respect to z. 



However, while the definition of as a total heat flux is intuitive and perfectly 
adequate for incompressible fluids (for which the theory was originally designed), a subtle 
but crucial problem emerges for compressible systems, where Tq^ 7^ 0. The total heat 
flux explicitly depends on Tq^, while the original system of equations ^ only knows about 
— Tq^, as discussed in Section 2.1 This remark suggests that the dynamically relevant 
quantity is instead 



kt\Toz - 



ptot ^ r.,,.,z J _ - ^ ^^^^ ^^2) 



ad I 

ad 



The flux thus defined, however, is no longer the total heat flux except when Tq^ = 0. 
Simplifying the resulting expressions for and -F*°* yields 

= 1 - T, + {wf) , 

= 1 - /ij + (^/i) , (13) 



These expressions are the ones actually used by Rosenblum et al. (2011). The system of 
equations ( 11 ) and ( 13 ) are now mathematically consistenlj^mean-field versions of the original 
system doj). 



We now define two non-dimensional quantities: 

1-T. 



•^An alternative, but equivalent, way to resolve the problem discussed here is to introduce the "potential 
temperature" i? commonly used in the atmospheric literature (e.g. Holton (1992)). The evolution equation 
for d is identical to that for T, except that the adiabatic gradient of -d is zero by construction. 
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^tot 

7tot = ■ (14) 

The first, Nu^, reduces to the much more commonly used temperature thermal Nusselt 
number (i.e. the ratio of the total heat flux to the diffused heat flux) when Tq^ = 0. In what 
follows, we call it the "thermal Nusselt number proxy". We also refer to the second, 7^"^, as 
the "flux ratio", for simplicity. When Tq^ = 0, it reduces to the total buoyancy flux ratio 
commonly used in physical oceanography. 



The theory then continues exactly as in Rosenblum et al. (2011), by assuming that Nut 



and 7t~J each depend only on the fluid parameters Pr and r and on the local inverse density 
ratio. The latter can vary with z as a result of the large-scale background temperature and 
compositional perturbations T and 

^-1 ^ /3(/io. + («//3)|To,-rof|/Zj ^ R^' -J, 
' a(To,-Tof + |To,-Tof|T,) 1 - ' 

where, for clarity, we first expressed R'^ as the ratio of dimensional quantities and then as 
the ratio of non-dimensional quantities. 



Combining (11), (14) and (15) yields a nonlinear system of equations describing the 



spatio-temporal evolution of the large-scale fields: 

dT _ (9F*°* 
dt dz 



i?0 - fJ'z 



where = NuT(i?;^ Pr, r)(l - T,) , i?^'^ = . (16) 

If 7t-t'(i?;^Pr,r) and NuT(i?;^Pr, r) are known, finite, non-zero, and smooth enough, 
then the system of equations is closed and well-posed. It has a trivial steady-state solution 
when Tz and JI^ are constant. This solution corresponds to a homogeneously, diffusively 
convective state with constant density ratio Rp^- Without loss of generality, we can choose 
our reference state T02, Tq^ and /ioz to be that steady-state solution, in which case R~^ = Rq^, 
and T2 = /i^ = 0. The flux ratio and thermal Nusselt number proxy of the homogeneous 
background turbulent state are noted as Nuq = Nut(-Ro^) %^ — 7tot(-Ro 



Solving (16) in the general case is numerically possible if the functions Nut and 7^"^ are 
known, but not particularly informative. However, we can linearize the mean-field equations 
around the previously defined homogeneously convective state, assuming that the large-scale 
perturbations T and JI have small amplitudes. To linear order, the local inverse density ratio 
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becomes 



(17) 



Noting that Nu^ depends on z via Rp^, it can be shown using the chain rule that, to hnear 
order, the temperature equation becomes 



dT 

~dt 



where 



A, 



^1 rfNu2 



dR-^ 



(19) 



while the linearized composition equation is similarly derived to be 



where 



-R, 



-1 dhtol) 



dR-^ 



(20) 



(21) 



Assuming normal modes of the form ~ g«K^+At^ finally get 

+ Ae [A2{1 - Ro%^) + Nuo(l - AiRo)] - fc^AiNugi?o = 



(22) 



This quadratic recovers the one obtained by Radko (2003) and Rosenblum et al. (2011) 



exactly, the only difference being in the physical interpretation of the quantities Nut and 
7t~t, as discussed above, when Vad 7^ 0. 



As originally discussed by Radko (2003), inspection of (22) shows that the condition 



for the existence of growing solutions is that the constant term in the quadratic should be 
negative, which only occurs when Ai is positive, i.e., when 'j^^l is a decreasing function of 
Rp^- In the diffusive case studied here, one can prove by inspection of the sign of the linear 



term in (22) that this sufficient condition is also a necessary condition for instability (by 
showing that even if there are complex conjugate roots to this equation, their real parts 



are negative). Radko (2003) also showed that the 7— instability theory suffers from an 
ultraviolet catastrophe whereby the mode growth rate is proportional to fc^ (so that modes 
with the smallest wavelengths always grow most rapidly). The theory, however, must break 
down when the layering mode wavelength becomes comparable with the basic instability 
wavelength. As a result, the actual mode which ends up growing out of the homogeneous 
turbulence is the one with the smallest wavelength for which the mean-field theory is still 
valid. Empirically, we find that the latter typically has a vertical wavelength that is about 
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2-4 times larger than the horizontal wavelength of the fastest growing mode of the basic 
instability according to linear theory (see Appendix A). In other words, the staircase typically 
forms with an initial step separation of about 25-50d. 

Finally, in order to identify more quantitatively the conditions for instability and predict 
its growth rate, we must measure the turbulent fluxes in the homogeneous phase of diffusive 
convection to estimate Jq^, Nuq, Ai and A2, for various values of Rq^, Pr and r. In all that 
follows, we therefore limit our definitions of 7^"^; and Nuy to the case where T = /I = 0, and 
so 

Nut = 1 + (wf) , 

7.;; - '-f^ . (23) 

1 + {wT) 

We also define for convenience a compositional Nusselt number 

Nu, = l+<||, (24) 

which measures the ratio of the total compositional flux to the diffused compositional flux 
in the homogeneous phase. With this definition, 

7.;:=rfl„-'^. (25) 



4. Measurements of the flux ratio 

As we have just shown in Section |3| double-diffusive layering is expected to occur spon- 



taneously whenever the flux ratio 7^"^ defined in (23) is a decreasing function of the inverse 
density ratio R^^ = V^/(V — Vad)- In what follows, we refer to the function ■J^oti^o^) 
as "the 7— curve". In order to establish when the 7— curve decreases and estimate the 7- 
instability growth rate, we now perform a series of numerical experiments, decreasing Pr and 
r down progressively towards the astrophysically relevant parameter regime, and measure 
both Nut(-Ro ^; Pr, r) 7tot (-^(7^! ''") whole range of density ratios unstable to 



diffusive convection (see ^4^). Section 4.1 describes our experimental setup and the manner 



in which we extract the flux ratio and the Nusselt numbers from the simulations. Section 



4.2 presents and discusses our results. 
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4.1. Experimental setup 



As in Rosenblum et al. (2011) and Traxler et al. (2011a), we use a computational box 
of size Lx = Ly = = lOOd, which is about 4-6 times the wavelength of the fastest- 
growing mode of instability, regardless of the parameters selected (see Appendix A). This 
domain size was found to be sufficiently large to yield statistically meaningful measurements 
of the turbulent fluxes while remaining computationally tractable in the increasingly extreme 
parameter regimes studied. 

We consider four values of Pr and r, equal to 0.3, 0.1, 0.03 and 0.01 respectively. The 
two smallest values, 0.03 and 0.01, are within the planetary parameter range. For each (Pr,r) 
pair, we run a number of simulations varying -Rg ^, selecting preferentially values close to one 
to capture the expected decreasing part of the 7-curve. Since the code is a Direct Numerical 
Simulation with no sub-grid model, each numerical experiment has to be fully resolved on all 
scales. Prior to each full-scale run, we test various resolutions and select the most appropriate 
one based on inspection of the vorticity, velocity and chemical composition field profiles and 
spectra. Runs with R^^ close to unity require the highest spatial resolution while runs with 
Rq^ close to marginal stability require lower spatial resolution, but much higher temporal 
resolution and longer integration times to follow simultaneously the buoyancy frequency 
timescale and the much slower instability growth and saturation timescales. Tables 1 and 2 
summarize the parameters selected and resolution for all our simulations. 

For each simulation, the PADDl code returns the non-dimensional instantaneous fiuxes 
integrated over the entire computational domain as diagnostics of the simulations. However, 
the thermal Nusselt number proxy and fiux ratio defined in the derivation of the 7— instability 
theory (see Section [3]) are only meaningful when viewed as temporal averages taken during 
a time where the system is in the assumed homogeneous, quasi-steady, diffusively convective 
state. Identifying that state, unfortunately, turns out to be significantly more difficult than 
expected. Figure |3] shows that transport in diffusive convection is much more variable than 
in the related fingering regime, where the layering theory and the methods for extracting 
small-scale fiuxes were first derived (Traxler et al. 2011a). The underlying reason for this 



difference actually remains to be determined. Our selected domain size, for example, was 



initially chosen by analogy with studies of transport in fingering convection by Traxler et al. 
(2011a) and Traxler et al. (2011b), where it was found to be "[...] small enough to suppress 



any secondary large-scale instabilities" (Traxler et al. 2011a). We find here, by contrast 



that large-scale perturbations (layers, large-scale gravity waves) do in fact grow even in such 
a small domain, and cause the observed variability in the fiuxes. 

In Appendix B, we discuss the problem in detail, and propose a systematic method to 
identify the homogeneous state described above, and extract the fiuxes, Nusselt numbers 
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Pr 


r 


D— 1 




^tot 


1 o 

layers : 


0.3 


0.3 


1.10 


384, 384 


907 


Y 


0.3 


0.3 


1.15 


384, 192 


1424 


Y 


0.3 


0.3 


1.20 


192,192 


1582 


Y 


0.3 


0.3 


1.25 


192, 192 


3251 


Y 


0.3 


0.3 


1.35 


192, 192 


2570 


? 


0.3 


0.3 


1.50 


96,96 


1999 


N 


0.3 


0.3 


1.60 


96,96 


1999 


N 


0.3 


0.3 


1.85 


96,96 


15000 


N 


0.1 


0.1 


1.10 


384, 384 


1114 


Y 


0.1 


0.1 


1.25 


384, 384 


1310 


Y 


0.1 


0.1 


1.50 


192, 192 


3095 


Y 


0.1 


0.1 


1.75 


192, 192 


3028 


? 


0.1 


0.1 


2.25 


192, 192 


3531 


N 


0.1 


0.1 


3.25 


192, 192 


4138 


N 


0.1 


0.1 


4.25 


192, 192 


10870 


N 


0.1 


0.1 


5.00 


192, 192 


22146 


N 


0.03 


0.03 


1.50 


576, 768 


2104 


Y 


0.03 


0.03 


2.00 


576, 576 


1587 


? 


0.03 


0.03 


2.50 


576, 576 


1311 


? 


0.03 


0.03 


3.00 


576, 576 


2215 


N 


0.03 


0.03 


4.00 


384, 384 


3148 


N 


0.03 


0.03 


5.00 


288, 288 


5929 


N 


0.03 


0.03 


10.00 


192, 192 


14845 


N 


0.01 


0.01 


1.50 


576, 576 


2987 


Y 


0.01 


0.01 


2.00 


576, 576 


4163 


? 


0.01 


0.01 


2.50 


576, 576 


1745 


? 


0.01 


0.01 


3.00 


576, 576 


2114 


N 


0.01 


0.01 


4.00 


384, 384 


2911 


N 


0.01 


0.01 


10.00 


288, 288 


8138 


N 



Table 1: Presentation of the various runs performed. The first three columns present the 
system parameters. All runs are in cubic domain of size (lOOrf)^. The resolution (in terms 
of equivalent mesh-points Nr^^y, N^) is always the same for the two horizontal direction, but 
occasionally differs in the vertical direction for runs that are expected to transition into 
layers. The total integration time is given in non-dimensional units as ttot- Finally, we 
indicate whether we see layers emerge or not. Runs with a question mark are runs for which 
we might expect layer formation based on the 7— instability criterion, and the actual position 
of the minimum of the curve, but where we have not seen evidence for it (yet). 
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Pr 


r 


i?0 ^ 




^tot 
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0.1 


1.20 
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0.1 


1.40 
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Y 


0.3 


0.1 


1.70 
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? 


0.3 


0.1 


2.00 
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N 


0.3 


0.1 


3.00 


192,192 


5506 


N 


0.1 


0.3 


1.10 


192, 192 


1759 


Y 


0.1 


0.3 


1.20 


240, 240 


1923 


Y 


0.1 


0.3 


1.30 


192, 192 


1702 


? 


0.1 


0.3 


1.50 


192, 192 


1946 


N 


0.1 


0.3 


2.00 


192, 192 


4314 


N 


0.3 


0.03 


1.10 


576, 576 


430 


Y 


0.3 


0.03 


1.25 


384, 384 


628 


Y 


0.3 


0.03 


1.50 


384, 384 


1052 


Y 


0.3 


0.03 


2.00 


288, 288 


937 


? 


0.3 


0.03 


3.00 


192, 192 


6847 


N 


0.03 


0.30 


1.10 


576, 576 


1574 


Y 


0.03 


0.30 


1.20 


384, 384 


2262 


Y 


0.03 


0.30 


1.35 


384, 384 


4100 


? 


0.03 


0.30 


1.50 


384, 384 


3177 


N 


0.03 


0.30 


2.00 


288, 288 


5750 


N 



Table 2: (Continued from Table 1.) 
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and flux ratio in that phase. The results are presented below. 



4.2. Nusselt numbers and flux ratio 

Our measurements for Nut and and 7^^, obtained using the method described in 
Appendix B, are summarized in Figures |4] and [s] respectively for each parameter pair (Pr, r). 
The full dataset is presented in Appendix B. 




100 



0.01 




0.01 



Fig. 4. — (a) The thermal Nusselt number proxy as a function of the reduced stability parameter 
r as defined in the main text, (b) The compositional Nusselt number as a function of the reduced 
stability parameter r. 

Figures |4^ and |4|d show Nuy — 1 and Nu^ — 1 respectively. Each curve represents one 
parameter pair (Pr, r), and is plotted against the stratification parameter r, where 



Rq ^ 



R-^ - 1 



(26) 



This quantity is introduced, following Traxler et al. (2011b), to re-map the instability range 
into the interval [0, 1], with r = corresponding to Ledoux criterion (r < being unstable 
to overturning convection) and r = 1 corresponding to the marginal stability limit (r > 1 
being fully stable). This new variable eases the comparison between the various datasets, 
and can be interpreted as a rescaled bifurcation parameter which measures the distance to 
stability/overturning instability. 



Figure 4 is reminiscent of a similar figure obtained by Traxler et al. (2011b) in the 



fingering regime. The thermal Nusselt number proxy is of the order of a few tens, and the 
compositional Nusselt number is of the order of a few hundreds for systems which are nearly 
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Ledoux-unstable. Both rapidly drop to one close to the marginal stability limit (r — )■ 1, 
Rq" —7- R^^\ Since a real Nusselt number can also be viewed (in the Boussinesq limit) as 
the ratio of the effective diffusivity (turbulent + microscopic) to the microscopic diffusivity, 
with 

Deff = Nu^K^, (27) 

our results show that turbulent compositional transport can be significant for more unstable 
systems. An equivalent interpretation for heat transport is more delicate, since Nut can 
only be viewed as a Nusselt number when Vad = 0. 




•^0 r 

Fig. 5. — (a, left) The flux ratio 7^~J obtained using the averaging methods discussed in Appendix 
B, for various values of Pr and r, as a function of ■ Only the interval R^^ G [1, 3] is shown to 
emphasize the region of decreasing 7^J: . The larger symbols indicate which runs eventually lead to 
layer formation, (b, right) The same results plotted against the instability parameter r as defined 
in the main text. The value of 7^"^: for r = 1 is the ratio of the diffusive fluxes, 7^^ (^ = f ) = ^o^'^- 

Figures [5^ andjsjj show the flux ratio 7^"^ measured in the simulations, as a function of 
Rq^ and as a function of r respectively. Both figures reveal many interesting features. We find 
that for all (Pr, r) explored, there exists a region where 7^^ decreases with Rq^ (equivalently, 
with r), hence, where layer formation is possible according to the 7— instability theory. We 
can therefore immediately compare our theoretical expectations with the actual outcome of 
the simulations: Figure [5^ and[5]D show runs which lead to layer formation as larger symbols. 
For larger values of Pr and r (i.e. Pr, r equal to 0.3 or 0.1), we confirm that layers indeed 
form whenever 7^"^ is a decreasing function of -Rq ^ hence validating the adequacy of Radko's 
criterion. For lower Pr and r, computational constraints limit our ability to validate Radko's 
theory as systematically as in the higher Pr and r case. Indeed, the layering mode growth 
rate depends on the derivative of 7tot(-Ro^) (see Section [3|, so the emergence of layers can 
be delayed significantly in runs with values of Rq^ close to the minimum of the curve (see 
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for example Figure gfor Pr = r = 0.1, Rq^ = 1.5). Since simulations at lower values of Pr 
and T require considerable spatial resolution, we were not always able to integrate them for 
enough time to see the emergence of layers. We have seen them for very low values of R^^ 
where the 7— curve decreases most rapidly, and expect that they should appear for slightly 
higher values of Rq^ as well. 

The fact that layering is possible in diffusive convection at low Pr and r is in stark 



contrast with results from the fingering regime (Traxler et al. 2011b), where 7tot always 
seems to increase with Rq in the same limit. This rather remarkable difference in behavior 
is actually fairly easy to understand. Indeed, let us first look at the behavior of the 7— curve 
close to marginal stability. In the corresponding runs, turbulent transport becomes negligible 
(see Figure |4]), so the flux ratio is dominated by diffusive transport. Mathematically speaking, 

T Nu t(1 — t) 

Nu., Nu„ ~ 1 ^ 7- = ^ ~ ' = + r , (28) 

which explains the observed oblique asymptote up to the limiting diffusive value tR~^ at 
r = 1 (see Figure [5}d). In the diffusive regime considered here, tR~^ is always smaller 
than one. However a similar argument applies in the fingering regime and yields 7tot = 
= -^-2 ^ i_ This limit "pulls up" the end of the 7— curve to very large values, 
effectively preventing the existence of a region where 7^"?; decreases with r. 



5. Theoretical predictions for the flux ratio 

The numerical simulations we have been able to perform sample parameter space rea- 
sonably comprehensively for Pr and r between 0.01 and 0.3, in particular for values of Rq^ 
close to unity. We found that for all parameter pairs (Pr, r) studied, there exists a interval 
Rq^ G [1, -Rl^] where the function 7tot (-^o ^) decreases, and that spontaneous layer formation 
indeed occurs in that region as expected from the 7— instability theory. However, in order 
to create a model for diffusive (semi-) convection that can be used practically and efficiently 
in a planetary or stellar evolution code, it would be preferable to have an analytical or semi- 
analytical theory for the position R^^ of the minimum of the 7— curve, rather than having 
to rely on interpolations or extrapolations of the available dataset presented in Tables 6 and 
7. In this section, we propose such a model. 
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5.1. Theoretical model for 7^^^ 



While we are looking for a model of the flux ratio 'Jt^u it is interesting to note that a 
method for estimating the turbulent flux ratio 



7turb 



{wf) 



(29) 



from linear theory was flrst proposed by Schmitt (|1979|) in the context of flngering convec- 
tion. 



Schmitt's theory adequately captures the shape of the curve 7turb('^) measured from 
laboratory ( Schmitt]|1979 ) and numerical experiments ( Traxler et al.|[2011b ), and in particu- 
lar its dependence on Pr and r, although the exact value of 7turb for a given value of r could 
be off by 20-40%. As such, it should be considered as a qualitatively accurate indicator of 
scalings and trends, but is quantitatively reliable only within factors of "a few" . 

Schmitt's method can straightforwardly be applied to diffusive convection, as derived 



in Appendix A3. The resulting expression for 7^^^^^^ is given by equation (A13), and depends 
on the growth rate and wavenumber of the most rapidly growing mode according to linear 
theory at the selected parameters Rq^, Pr and r. The latter can be found numerically quite 
easily by solving simultaneously a cubic and a quadratic equation. If is known, then 



-1 _ T-Rq + 7tm.b 
/tot 



-1 



[wT) 



^^o' + 7tuib(NuT 



(wT) 



(Nut 



(30) 



where we have used equation (23) to express the turbulent heat flux in terms of Nur- AH 



that remains to do is to create a model for Nu-r as a function of the system parameters Rq ^, 
Pr and r. 



We now return to the results of the numerical simulations presented in Section |4.2[ We 
find that we can satisfactorily fit the behavior of Nut — 1 for large values of r provided 
Nut — 1 oc (1 — r). Close to overturning instability, on the other hand, we find that a first 
satisfactory fit to the data has Nut — 1 oc (1 — t)/{Rq^ — 1). This functional dependence 
is not unexpected, since Rq^ — 1 is the non-dimensional background density gradient, and 
since diffusive convection relies on r 7^ 1 to operate, and is much more efficient the smaller 
the value of r. Combining this fit with the large Rq^ limit suggests a functional form with 
Nut-1 oc (l-r)(l-r)/(i?o^-l)- One can fit the proportionality constant for runs with Pr 
= r, and obtain a rather good match to the data. However, the resulting expression is less 
satisfactory for Pr 7^ r. Further investigation reveals that an even better fit can be obtained 
with 



{wf) = Nut - 1 = (0.75 ± 0.05) 



pr\0.25±0.15 



R, 



1-r) 



(31) 
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The large uncertainty on the power index of the term (Pr/r) comes from the uncertainty on 
the measurements themselves, compounded with the short range of Pr/r values available. 
However, its exact value does not matter much for the 7^^^^ predictions. 
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Fig. 6. — Comparison between expression (31) and the data presented in Tables 6 and 7 



Figure ^ compares our empirical fit for Nut — 1 given by equation (31) to the actual 



data. This fit is satisfactory for our current purposes, although we recognize that a better 
theoretically-motivated one should be sought in the future if we wish to improve on the 
model further. 



Using (A13), (30) and (31) we can now estimate 7t~J semi-analytically. Figure [t] com- 
pares our predictions with the data presented in Figure [5]d. As expected from the limitations 
of Schmitt's method, and uncertainties in our fit for the turbulent heat flux, the model does 
not match the data perfectly. Generally speaking, we find that the predicted value of r at the 
minimum is somewhat overestimated by the model, by 20-40%. The slope of the 7— curve is 
thus also affected. These discrepancies are larger and/or more apparent for runs with larger 
values of the Prandtl number, for which the position of the minimum occurs for larger values 
of r. However, it is nevertheless rather remarkable to see how well our model accounts for 
the shape of the 7— curve, and in particular the variation of the position of the minimum 
with Pr and r. The value of 7t~J at the minimum is also robustly predicted by the model. 
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r 

Fig. 7. — Comparison between tlie model prediction for 7tot(r) and the data presented in Tables 6 
and 7 and Figure [5] The same color scheme is used in this figure to represent the various parameter 
pairs (Pr, r) as in Figures 6, 7, and 9. 

5.2. Model trends 

Using the model described in the previous section, we can now estimate the value of the 
inverse density ratio R~[^ for which 'j^^l is minimal, for a range of parameter values beyond 
those for which we were able to run numerical simulations. The results are presented in 
Figure [sj for Pr and r varying from 10"'' to 1, and show contour plots of Rj^^ (bottom) and 
of the corresponding = {Rj} — — 1) (top). Note that both tl and estimated 

directly from the model, as shown here, are likely to overestimate the true position of the 
minimum by about 20%-40% (see previous section). 

Overall we find that the relative fraction of the total instability range unstable to layering 
decreases as Pr and r decrease (i.e. the value of r at the minimum of the 7— curve decreases). 
However, since the instability range itself increases as Pr and r tend to zero, the value Rl^ 
below which layers can spontaneously emerge actually increases significantly. We find that 
it is of the order of a few for planetary values of the diffusivities, and of the order of a 
few hundreds to a thousand for the stellar parameter regime. The layering instability, and 
its implications on increasing the heat and compositional transport properties of diffusive 
convection, is thus likely to play an important role in stellar and planetary astrophysics. 
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Fig. 8. — Contour plots of the position of the minimum of the 7— curve as a function of Pr and 
T, written as R~[^ (bottom) and ri = {R^} — 1)/{R^'^ — 1) (top). Note that both and Rj} 
calculated directly from the model, as shown here, are likely to overestimate the true position of 
the minimum by about 20%-40%. 

6. Comparison of the layer growth rates with theory 



Rosenblum et al. (2011 ) already showed that Radko's 7— instability theory, when applied 



to the case of diffusive convection, correctly accounts for the growth rate of the emergent 
staircase in their simulations at Pr = r = 0.3 with R^^ = 1.2. In this section, we check 
that this is still true at lower values of Pr and r, and compare two methods for estimating 
the mode growth rate: one based on the experimentally-determined functions 7^^ (Rq^) and 
Nut(-Ro ^) listed in Tables 6 and 7, and one based on the model functions proposed in Section 

m 

As described in Section |3] the growth rate A of a layering-mode with vertical wave- 



number k is the solution of the quadratic (22). Estimating A thus requires first estimating 



Nuo = Nut(-Ro^) 7o ^ = 7tot(-Ro^) respectively, as well as the derivative terms Ai and 



A2 defined in equations (21) and (19) respectively. This can either be done using the actual 



experimental data, or using our new semi-analytical theory (see Section [s]). When using the 
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experimental data, Ai and A2 are calculated using either one-sided or two-sided derivatives, 
depending on the data points available. When using the semi-analytical model functions, Ai 
and A2 are always calculated using two points on either sides of the selected value of Rq^. 

To illustrate the process, we compare the 7— instability theory with the data from the 
Pr = r = 0.03, Rq^ = 1.5 run. Table [s] shows the results of our estimates for the layering 
mode growth rate A using the two different methods. The experimentally-derived results 
are expected to be more accurate since they do not rely on any modeling. Reassuringly, 
however, we find that the model-derived growth rate is within 20% of the experimentally 
derived one. It is interesting to note that while the model-estimate for A2 seems to be off by 
an order unity, this discrepancy does not affect the growth rate estimate much. This remark 
is valid for many of the cases studied (although in many other cases A2 is well-predicted by 
the theory). 

Experiment (Ml) Model (M2) 



Nuo 2.36 2.41 

7o"^ 0.31 0.38 

Ai 0.33 0.49 

A2 2.34 4.36 

A/k^ 0.31 0.36 



Table 3: Layering mode growth rate using Nuq, 7o ^, Ai and A2 from the experimental 
data (Ml) and from the model presented in Section [s] (M2) respectively, for the run with 
Pr = r = 0.03, R^^ = 1.5. Note that since A is proportional to k^, we list the proportionality 
constant A/k"^ for more generality. 



Figure [9] compares our estimates for A with the actual mode growth observed in the 



simulations. As shown by Stellmach et al. (2011) and Rosenblum et al. (2011), a convenient 



way of extracting the amplitude of the layering mode is to look at the Fourier expansion 
of the density field, and isolate the mode with zero horizontal wavenumber, and a vertical 
wavenumber k = n{2'K/Lz) where n is the number of steps in the emergent staircase. Figure 
[9] shows the square of the norm of that mode (i.e. its power IpnP), and compares it with an 
exponential function proportional to e^^* (the normalization being arbitrary). Both growth 
rate estimates correctly account for the observed mode growth, with the experimentally- 
derived growth rate faring somewhat better, as expected. However, it is reassuring to see 
that the model proposed in Section |5] works quite well too. 

Once the mode's amplitude grows beyond a certain critical value, its density profile is 
no longer monotonously decreasing. When this happens, localized regions become unstable 
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to overturning convection, and a fully-formed staircase rapidly appears. The threshold for 



overturning instability for a mode with n steps was calculated by Rosenblum et al. (2011) 
to be, in terms of its density spectral power. 



\P 



2 

n I conv 



1 - ^0 ^ 
2n 



(32) 



Figure [9] clearly shows that the mode growth rapidly stops after its amplitude crosses that 
threshold. 



Pr=0.03, T=0.3, Rq =1.5 
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Fig. 9. — Comparison between the model prediction for the layering mode growth rate and the 
actual data, for Pr = r = 0.03, R^^ = 1.5. The emergent mode observed has n = 3 steps so 
k = 3{2tt/Lz)- The case using the growth rate derived from the experimental data is shown as Ml, 
while the one using the growth rate derived from the model for the thermal Nusselt number proxy 
and flux ratio proposed in Section [5] is shown as M2. Also shown is the critical amplitude for onset 
of overturning convection, as a horizontal line. The mode growth notably changes upon reaching 
this amplitude, and quickly saturates after that. 



Applying the same method to all runs which eventually result in layer formation, we 
find that the layer growth rate predicted by the solution of equation (22) always correctly 
accounts for the observed mode growth in the simulations. Furthermore, while the growth 
rate predicted from experimentally-derived values of Nuq, 7(7^, and A2 is always better 
than the model-derived ones, the latter are nevertheless satisfactory estimates too, and are 
always within 10-30% of the correct value. These results complete the validation of Radko's 
theory, as well as our model estimates for 7^"^ and Nut- 
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7. Conclusion and prospects 

The ultimate goal of this series of papers is to propose a new model for mixing by 
diffusive convection. In the work presented here, we ran and analyzed a very extensive 
suite of numerical simulations of the process in a wide range of parameter space. We have 
shown that in the astrophysically-relevant low Prandtl number (Pr = z//kt), low diffusivity 
ratio (r = n^/ kt) regime, diffusive convection can take one of two forms depending on the 
local inverse density ratio Rq^ = V^/(V — Vad): moderately efficient layered convection 
at lower Rq^, or inefficient wave-dominated "oscillatory" convection for higher Rq^ (see 
Figure [T|. We have confirmed through numerical and analytical work that a spontaneous 
transition into layered convection occurs, under predictable circumstances, through a linear 
mean-field instability of the initial state of oscillatory convection. This instability is called 
the 7— instability. It was originally suggested in the oceanographic context of fingering 



convection by Radko (2003) and later applied to double-diffusive convection in astrophysics 



by Rosenblum et al. (2011) 



Since transport in layered convection is much more efficient than in the non-layered 
case, a crucial element of any new model of diffusive convection will be the availability of a 
practical criterion for determining, for given Rq^, Pr and r, whether a system is expected 
to transition into layers or not, and on what timescale. The original 7— instability theory 
provides such a criterion, but the latter can only be used in practice provided experimental 
measurements of the turbulent fluxes in homogeneous diffusive convection, for the same 
parameters, are available. We provide such measurements here for the planetary parameter 
regime, but similar results are unlikely to ever be available for the much more extreme stellar 
parameter regime. 

Based on these considerations, we then proposed a new empirically motivated model for 
the turbulent fluxes, which enable us to derive a completely parameter-free semi-analytical 
criterion to determine, for any given fluid in the astrophysical regime (Pr, r <^ 1) and 
any given stratification (R^^), (a) whether a system is expected to transition into layers 
or not, and (b) on what timescales the layers are expected to emerge. Our model was 
found to fit the available numerical data very well, and can therefore be used very reliably 
within the same region of parameter space (e.g. the planetary parameter regime, for which 
Pr, r ~ 10^^ — 10^^). We further propose that it should also be used in regions of parameter 
space for which fully resolved simulations are not available, namely in the stellar parameter 
regime (Pr, r ~ 10"^ — 10~^). Based on this model, we find that layered convection is 
theoretically expected in stellar interiors for a fairly wide range of parameter space, with 
Rq^ between 1 and about 1000. 

Our results answer, at least approximately, the first of the three questions we initially 
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raised: (1) Under which conditions do layers form? (2) What is the transport rate in layered 
convection and (3) What is the transport rate in non-layered "oscillatory" convection. In 
subsequent papers in this series we will continue our investigation by answering questions 
(2) and (3). 

We thank Nic Brummell, Jonathan Fortney, and Douglas Gough for fruitful discussions. 
G.M., P.G. and T.W. were supported by funding from the NSF (NSF-0807672), and benefited 
from the hospitality of the ISIMA program during the summer of 2011. A.T. P.G. and T.W. 
were funded by the NSF (NSF-0933759). Part of the computations were performed on the 
UCSC Pleiades supercomputer, purchased with an NSF-MRI grant. Others used computer 
resources at the National Energy Research Scientific Computing Center (NERSC), which 
is supported by the Office of Science of the US Department of Energy under contract DE- 
AC03-76SF00098. Figure 1 was rendered using ViSiT. P.G. thanks LLBL Hank Childs for 
his excellent support of the software. 

A. Appendix: Linear stability of semi-convection, asymptotic solutions for 

low Pr and r, fastest-growing modes 

The system of equations ^ can be linearized and solved for the fastest growing modes of 
diffusive convection. These linear solutions can then be studied further to obtain asymptotic 
scalings at very low Pr and r, and to derive predictions for the turbulent buoyancy flux ratio 
(see Section [sj. 

A.l. Linearized equations for the fastest-growing modes 

We first linearize ^ around T = p, = and li = 0, assuming all perturbations are 
normal modes of the form q = ge«'2;+imy+«fc2+At -^j-^g^-g g ^ {T, /i,u}. Hatted quantities are 
now the amplitudes of the perturbations, while I and m are the horizontal wave-numbers, k 
is the vertical one, and A is the growth rate. The latter are all non-dimensional. 

We are interested in the fastest-growing modes only, which can be shown to have k = 

. They correspond to purely vertical 
fluid motions. They are rotationally invariant around the vertical direction, so without loss 
of generality we can align the horizontal wavenumber with the x-axis choosing m = 0. After 



as in the flngering case (|Radko|2003j|Traxler et al.|2011b| 
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some simplifications, the resulting system of equations for the mode amplitudes are 



XT-w = -I'T , 

\w = -Fil^w + Pr(f - il) . (Al) 

It has a non-trivial solution provided the growth rate A satisfies the following cubic equation: 
A 



Pr 



+ (A + /^)(A + Tl') - (A + rl') + R,\X + l') = . 



(A2) 



In the regime of interest, this cubic has one negative real root and two complex conjugate 



roots (Baines & Gill 1969). It can easily be shown that the complex conjugate roots (A = 

(A3) 



\r + iXi with A/ 7^ 0) satisfy 

A^ = 3A^ + 2P\r (r + Pr + 1) + (r + Prr + Pr) + Pt{R^' - 1) , 
and that A^ satisfies the cubic 

8A^ + SPxUt + Pr + 1) 
+2\n (r + Prr + Pr + (r + Pr + 1)') + Pr(i?o ^ - 1)] 
+1%T + Pr)(r + l)(Pr + 1) + /^Pr(i?o ^ (r + Pr) - (Pr + 1)) = . 

The latter has a positive solution if and only if Rq^ E [1, where R~^ = 

The fastest growing modes are determined by fixing Pr, r and Rq^ within the instability 



(A4) 



range, and finding the value of / for which A/j is maximum by solving (A4) in conjunction 
with ^f- = 0, or in other words 

8A|(r + Pr + 1) + AXrP (r + Prr + Pr + (r + Pr + 1)') 
+3/^(r + Pr)(r + l)(Pr + 1) + Pr(i?o ^ (r + Pr) - (Pr + 1)) = . (A5) 

In what follows, we study the behavior of the solutions as a function of the reduced 
stratification parameter r, defined in ([26|). Note that, with this new variable, we have 

(A6) 



i?o ' (r + Pr) - (Pr + 1) = (r - 1)(1 - r) . 



A. 2. Asymptotic solutions at low Pr and r 



In general, one needs to solve (A4) and (A5) numerically to find the fastest growing 



modes for given Pr,r and Rq^- Here, however, we are interested in deriving asymptotic 
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solutions for low Pr and low r, since this is the parameter regime relevant for planetary and 
stellar interiors. In particular, we want to study how the growth rate and the wavenumber 
of the fastest-growing modes scale with these governing parameters. 



The solutions to equations (A4) and ^A5v can easily be found numerically, and the results 



are shown in Figure 10 for decreasing Pr (here with Pr= r). We see that the real part of the 
fastest growing mode's non-dimensional growth rate, Amax('^), appears to be proportional to 
Pr, while the corresponding horizontal wavenumber, lmax{f), becomes independent of Pr as 
Pr decreases. 

This behavior suggests a new rescaling of the governing equations to capture the asymp- 
totic regime (Pr, r — )■ 0) of the instability: Amax = PrA where A ~ 0(1), and Zmax = i where 
/ ~ 0(1). We also define = T/Pr, and assume that is order unity. Using (26) and keeping 
only the lowest terms in Pr, equations (A4) and (A5) reduce to the simple universaj^ system: 



2 + 



+ 1 



4PA + 3/^(0 +l) + (r-l) 



A + r(0+l) + r(r-ll 







(A7) 



This system still needs to be solved numerically for A and /, but only once for each value 
of and r. Figure 10 compares the solution of (A7) with the ones obtained by direct 
numerical solution of (A4) and (A5) for Pr = 10"*^ and = 1, and confirms our numerical 



and semi-analytical results. 

This linear asymptotic analysis helps us estimate the growth rate of this kind of double- 
diffusive instability for a broad range of parameters and determine the size of the basic 
unstable structures we are interested in. Dimensionally speaking, our results imply that the 
true lengthscale of the instability should always be of the order of a few d, where d was defined 
in and the growth rate of the instability should be of the order of PrKr/c?^ = VPtN 
where N is the thermal buoyancy (Briint-Vaisala) frequency. This information is useful for 
two reasons: first, to quantify the expected lengthscales or timescales in the real systems (i.e. 
stellar and planetary interiors), and secondly, to get some insight into the correct domain 
size and timestep to use in the numerical simulations. 



■^Note that this asymptotic Hmit is not uniformly vaUd for r — >■ 0. 
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A. 3. Semi-analytical prediction for the turbulent buoyancy flux ratio 

Let us consider the turbulent flux ratio 

7turb - 7T;f7 l^^J 
{wT) 



where (■) denotes a spatial average over the entire computational domain. Schmitt (1979) 
showed that it is possible to estimate this quantity for fingering convection using the velocity 
field, temperature and chemical composition perturbations corresponding to the linearly 
fastest-growing mode of instability. Since the unknown amplitude of the perturbations in 
this turbulent ratio cancels out, the remaining expression only depends on the known shape 
and growth rate of the perturbations. Here, we apply the same technique to estimate the 
turbulent flux ratio in diffusive convection. 



From the system of equations (Al), we see that the amplitudes of the vertical velocity. 



temperature and compositional perturbations of a given mode are related via 



A + /2 ' 

a = rA9) 

In order to calculate the fluxes, we must remember that, in the process of the linear analysis, 
the various flelds w, T and jl were deflned as complex variables, e.g. from q = qQ^^^+'^^'V+i'^^+^t- ^ 
under the implicit understanding that only their real parts are physically meaningful. Hence, 
in practice, 

_ mw)^m 

^''^^ (3?(w)3?(f)) ■ ^ ' 

Without loss of generality, w can be selected to be real so that 

3f?(w) = we^«* cos(/x + Xit) . (All) 



Then, using (A9), we flnd that 



^(^) =7^ -—^:—-^[cos{lx + Xit){\R + P) + sm{lx + \it)\i] , 

[■^R + I- ) + -^l 

= IX ,° n^2 , ^2 [cos(/x + \it){\R + rf) + sin(/x + A/t)A/] . (A12) 



Finally, forming the turbulent flux ratio and integrating the relevant quantities over the 
computational domain and over short timescales (i.e. over at least one oscillation period of 
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the basic instability), we get 7turb(^' ^ given mode with wavenumber / and growth rate 

\ = Xr + iXj as 



{Xn + rPy + Xj Xn + P 



(A13) 



where Xr and A/ are related via (A3). Applying this formula to the most rapidly growing 
mode, with wavenumber /max and growth rate Amax calculated in Appendix Al, yields the 
required estimate for the inverse turbulent flux ratio in our simulations. 



B. Appendix B: Extraction of mean fluxes and results 

B.l. Protocol for extracting mean fluxes and measuring 7^"^ from the 

simulations 

In what follows, we describe our protocol for measuring the mean turbulent fluxes in the 
homogeneous phase of diffusive convection (prior to the emergence of large-scale structures). 
This involves first creating a systematic method to identify the start and end times [tstart, ^end] 
of this phase and then estimating the fluxes and related errorbars. 



B.1.1. Selection of tstart- 

As seen in Figure [3} the turbulent flux typically peaks then drops quite sharply during 
the saturation of the primary instability, and then grows more slowly towards its value 



in the homogeneous double-diffusively convecting state. As shown in Figure [TT| the same 
description applies to the behavior of the total kinetic energy in the system. It thus appears 
that the system needs a little bit of time to "recover" from the saturation. In order to extract 
meaningful averages, we therefore need to select the start of the averaging process well-past 
the main saturation peak. We also need to define tstart in a manner that is meaningful across 



all simulations. Figure 11 illustrates our process: we define first the "width" of the saturation 
peak At as illustrated, and then choose tstart accordingly, about 2At past the peak. While 
this choice is arguably somewhat arbitrary, it does satisfy the requirements listed above. 
Furthermore, the estimated values of 7^"^ are not particularly sensitive to the choice of tgtart 
as long as it is indeed well-past the saturation peak. 
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B.1.2. Selection o/tend in the non-layered case 



As noted earlier, by contrast with Rosenblum et al. (2011) we find that even in the 



non-layered case the system does not necessarily remain in a state of homogeneous, small- 
scale diffusive convection but sometimes becomes dominated by larger-scale coherent gravity 
wave ^ While the precise reason for the emergence and synchronization of these waves re- 
mains to be determined, their associated dynamics lead to a rather different type of transport 
than in more homogeneous diffusive convection. For this reason, we must identify when the 
waves first "take over" and restrict our measurements of the turbulent fiuxes prior to that 
time. 



Shown in Figure [12] is the total kinetic energy in the simulation, as well as the total 
kinetic energy in the six highest-amplitude families of gravity wave modes. By "families", 
we imply the following. A single gravity-wave mode, in this triply-periodic simulation, can 
be identified with the Fourier mode proportional to exp^ik^x + ikyi/ + ikzz) , where {k^, ky, k^) 
is the mode wave-vector. A "family" of modes is defined as the ensemble of all the modes 
with the same geometry given the symmetries of the system, i.e. the same values of \kz\ and 
the same values of \kh\ = ^Jk^. + k'^. In what follows, we classify the modes for simplicity of 
notation based on their periodicity: the single mode {0,2, —1} for example corresponds to 
one with k^ = 0, ky = 2(2'K/Ly), k^ = —{2tt/Lz). The family of modes 021 then corresponds 
to an ensemble of 8 modes: {0, 2, 1}, {0, 2, -1}, {0, -2, 1}, {0, -2, -1}, {2, 0, 1}, {2, - 1}, 
{—2, 0, 1} and finally {—2, 0, —1}. Finally, the total kinetic energy in the mode family is just 
the sum of that of the individual modes. 



Figure 12 shows that the evolution of the total kinetic energy of the system is very 
similar to that of the turbulent fiuxes for the same simulation (see Figure |3]): an extended, 
apparently quasi-steady turbulent state between t ~ 600 and t ~ 2000, followed by a wave- 
dominated phase. It also reveals that the family of modes which dominates the system 
beyond t = 2000 is the 012 family, and that the strong oscillatory signal in the total kinetic 
energy (and the turbulent heat fiux) appears when the total kinetic energy in that single 
family exceeds half the total kinetic energy of the system (shown as the thin black line). 

We used a similar method to analyze every single simulation among the ones presented 
in Tables 1 and 2, comparing the total kinetic energy to that of various families of modes, 
and found that a robust (albeit empirical) criterion for determining the time tgw when a 
system becomes dominated by gravity waves is simply that the total kinetic energy in any 



Rosenblum et al. 



(2011 1 did not notice the emergence of the waves in their simulations, although a more 



careful re-analysis of their results shows that they were indeed present in some of the higher i?Q ^ runs. 
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given family of modes exceeds half the total kinetic energy in the system. In the non-layered 
case, we therefore take the "end-point" tend of the temporal average to be tend = tgw In 
some cases with Rq^ close to marginal stability, it can happen that the start and end times 
thus selected have tstart > ^nd- When this is the case, we discard the simulation (for the 
purpose of estimating the turbulent fluxes and their ratio). 



B.1.3. Selection of tend in the layered case 

Applying the method described in the previous section we find that, in runs which 
eventually show the emergence of a staircase, gravity-wave modes never dominate the system. 
However, since we are interested here in the process which leads to layer formation, extracting 
the flux ratio in the layered case is only meaningful prior to the formation of the first layers. 
So, whenever layers appear in the simulations, we set tend to be the time where the first set 
of layers appears. 



B.1.4- Averaging method and error estimates 

Once the relevant time interval has been determined, we need to measure the mean 
turbulent fluxes, construct Nut, Nu^ and 7t~J: and estimate our experimental error. For this 
purpose, we use a "4-intervals" method: we first divide the integration domain previously 
defined into four sub- intervals, and calculate the mean fluxes and therefore Nut, Nu^ and 



7t~t in each one of them according to (23) and (24). The final adopted value of Nut, Nu^ 
and 7t~t respectively is then the average of the four computed values, while the error is their 
standard deviation. The reason for using this method is clarified in the examples below. 

Let us first illustrate our procedure on the data from the simulation shown in Figure 
|2^, i.e. for the run that leads to layer formation (with Pr = r = 0.03, -Rg^^ = 1.5). We 
first estimate the start- and end-times of the homogeneous phase to be tstart = 775 and 
tend = 1200. The mean Nu^, Nu^ and 7^"^ in each sub-intervals are given in Table [ij as 
well as their final values and corresponding errorbars. These results illustrate the reason for 
using such a method to estimate the measurement "error" rather than a simple average over 
a single interval: the mean Nut and Nu^ increase steadily from one sub-interval to the other, 
showing that the system is not actually in a statistically quasi-steady state (as assumed by 
the 7— instability theory). However, this clearly does not prevent the layering modes from 
growing anyway, and as shown in Section [6} the theory still adequately accounts for their 
growth rate despite the non-stationarity of the homogeneous phase. As such, we have to do 
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the best with the data we have, and report on the values of Nu^, Nu^ and 7t~J accordingly, 
albeit with large errorbars which account for the slow temporal evolution of the system from 
saturation to the emergence of the staircase. 





^start 


^cnd 


Nut 


Nu^ 


7tot 


Interval 1 


775 


881.26 


2.082 


12.567 


0.272 


Interval 2 


881.26 


987.5 


2.256 


14.531 


0.290 


Interval 3 


987.5 


1093.75 


2.384 


16.519 


0.309 


Interval 4 


1093.77 


1200 


2.703 


20.448 


0.341 


Total 


775 


1200 


2.36 ±0.23 


15.9 ±2.9 


0.31 ±0.03 



Table 4: Illustration of our data averaging method for the run with Pr = r = 0.03, Rq^ = 1.5. 

Applying this method to the non-layered run shown in Figure 12 (with Pr=r = 0.1, 
Rq^ = 1.75) we find that the start and end of the homogeneous period are t = 650 and 
t = 2100, and the corresponding Nuy, Nu^ and 7^"^ computed are shown in Table [s] In this 
case the run is more stationary overall, leading to much smaller errorbars. 





^start 


^end 


Nut 


Nu^ 


7tot 


Interval 1 


650 


1012.5 


1.776 


3.279 


0.323 


Interval 2 


1012.5 


1375 


1.639 


2.830 


0.302 


Interval 3 


1375 


1737.5 


1.673 


2.992 


0.313 


Interval 4 


1737.5 


2100 


1.791 


3.308 


0.323 


Total 


650 


2100 


1.72 ±0.07 


3.10 ±0.20 


0.32 ±0.01 



Table 5: Illustration of our data averaging method for the run with Pr = r = 0.1, flg ^ = 1.75. 



B.2. Summary of the results 



The results of our analysis are summarized in Tables 6 and 7. 
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Pr 


T 




r 


start 


^end 


7tot 


Nut 




0.3 


0.3 


1.1 


0.09 


356 


506 


0.69 ±0.03 


8.43 ±2.02 


17.7 ±5.1 


0.3 


0.3 


1.15 


0.13 


370 


700 


0.61 ±0.01 


4.13 ±0.14 


7.26 ±0.38 


0.3 


0.3 


1.2 


0.17 


450 


920 


0.58 ±0.01 


3.21 ±0.20 


5.14 ±0.41 


0.3 


0.3 


1.25 


0.21 


450 


2000 


0.54 ±0.01 


2.50 ±0.21 


3.62 ±0.39 


0.3 


0.3 


1.35 


0.30 


550 


660 


0.51 ±0.01 


1.84 ±0.03 


2.33 ±0.07 


0.3 


0.3 


1.5 


0.43 


700 


950 


0.52 ±0.01 


1.53 ±0.04 


1.78 ±0.05 


0.3 


0.3 


1.6 


0.51 


940 


1070 


0.53 ±0.01 


1.37 ±0.02 


1.51 ±0.06 


0.3 


0.3 


1.85 


0.73 


2200 


2830 


0.57 ±0.01 


1.18 ±0.01 


1.21 ±0.01 


0.1 


0.1 


1.1 


0.02 


400 


480 


0.62 ±0.05 


8.92 ±2.31 


50.6 ±16.7 


0.1 


0.1 


1.25 


0.06 


500 


720 


0.47 ±0.02 


3.92 ±0.26 


14.7 ± 1.4 


0.1 


0.1 


1.5 


0.11 


575 


2150 


0.36 ±0.01 


2.21 ±0.10 


5.24 ±0.39 


0.1 


0.1 


1.75 


0.17 


650 


2100 


0.32 ±0.01 


1.72 ±0.07 


3.10 ±0.20 


0.1 


0.1 


2.25 


0.28 


820 


2700 


0.32 ±0.01 


1.43 ±0.05 


2.01 ±0.11 


0.1 


0.1 


3.25 


0.50 


1780 


2050 


0.36 ±0.01 


1.19 ±0.01 


1.32 ±0.02 


0.1 


0.1 


4.25 


0.72 


4300 


4900 


0.43 ±0.01 


1.05 ±0.01 


1.06 ±0.01 


0.03 


0.03 


1.5 


0.03 


775 


1200 


0.31 ±0.03 


2.36 ±0.23 


16.0 ±2.9 


0.03 


0.03 


2 


0.06 


1000 


1600 


0.20 ±0.01 


1.58 ±0.05 


5.30 ±0.52 


0.03 


0.03 


2.5 


0.09 


620 


1250 


0.19 ±0.01 


1.41 ±0.07 


3.53 ±0.35 


0.03 


0.03 


3 


0.12 


1300 


2215 


0.19 ±0.02 


1.35 ±0.07 


2.86 ±0.38 


0.03 


0.03 


4 


0.19 


1650 


2250 


0.19 ±0.01 


1.21 ±0.05 


1.88 ±0.19 


0.03 


0.03 


5 


0.25 


3400 


4800 


0.22 ±0.01 


1.23 ±0.05 


1.78 ±0.17 


0.03 


0.03 


10 


0.56 


5700 


6500 


0.30 ±0.01 


1.02 ±0.01 


1.03 ±0.01 


0.01 


0.01 


1.5 


0.01 


1230 


1840 


0.25 ±0.02 


1.95 ±0.14 


32.0 ±5.3 


0.01 


0.01 


2 


0.02 


1450 


3400 


0.19 ±0.01 


1.69 ±0.08 


15.9 ± 1.5 


0.01 


0.01 


2.5 


0.03 


1050 


1745 


0.13 ±0.01 


1.38 ±0.07 


7.39 ± 1.02 


0.01 


0.01 


3 


0.04 


900 


2200 


0.12 ±0.02 


1.31 ±0.09 


5.27 ± 1.11 


0.01 


0.01 


4 


0.06 


1150 


2911 


0.12 ±0.02 


1.27 ±0.08 


3.91 ±0.71 


0.01 


0.01 


10 


0.18 


3990 


5650 


0.13 ±0.01 


1.07 ±0.03 


1.35 ±0.14 



Table 6: Summary of the results. The first three columns are the run parameters, corre- 
sponding to those presented in Table 1. The fourth column shows the stability parameter 
r defined in (26). The 4th and 5th columns show the start and end times for the temporal 
average, as discussed in Section |4j The 6th, 7th and 8th columns show the fiux ratio 7^"^ 
and Nusselt numbers Nut and Nu^, as defined in equations (23) and (24). Three significant 
digits are shown for Nut and Nu^, and two for 7^"^. 
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ouu 


^ou 
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o.ou HZ u. lo 
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0.3 


0.1 


1.4 


0.18 


300 


400 


32 ± 02 


2 87 ± 06 


6 58 ± 52 


0.3 


0.1 


1.7 


0.31 


460 


1200 


0.26 ±0.01 


1.78 ±0.05 


2.78 ±0.05 


0.3 


0.1 


2 


0.44 


620 


920 


0.25 ±0.01 


1.42 ± 0.04 


1.81 ±0.12 


0.1 


0.3 


1.1 


0.06 


480 


620 


0.68 ±0.01 


4.52 ± 0.46 


9.29 ±1.10 


0.1 


0.3 


1.2 


0.11 


650 


1100 


0.61 ± 0.02 


2.80 ±0.21 


4.79 ± 0.49 


0.1 


0.3 


1.3 


0.17 


660 


1650 


0.56 ±0.01 


1.92 ±0.07 


2.76 ±0.14 


0.1 


0.3 


1.5 


0.29 


850 


1180 


0.57 ±0.01 


1.62 ±0.05 


2.05 ±0.08 


0.1 


0.3 


2 


0.57 


1500 


2350 


0.63 ±0.01 


1.17±0.01 


1.23 ±0.01 


0.3 


0.03 


1.1 


0.03 


225 


325 


0.57 ±0.02 


19.3 ± 1.9 


332 ± 42 


0.3 


0.03 


1.25 


0.09 


190 


500 


0.36 ±0.06 


6.62 ±2.26 


63.3 ±33.4 


0.3 


0.03 


1.5 


0.17 


220 


930 


0.20 ±0.02 


2.98 ±0.43 


13.0 ±3.1 


0.3 


0.03 


2 


0.34 


270 


937 


0.12 ±0.01 


1.64 ±0.10 


3.17 ±0.35 


0.03 


0.3 


1.1 


0.05 


900 


1420 


0.72 ± 0.01 


4.63 ± 0.32 


10.1 ±0.9 


0.03 


0.3 


1.2 


0.09 


1100 


1850 


0.68 ±0.05 


3.42 ± 1.36 


6.48 ±3.10 


0.03 


0.3 


1.35 


0.17 


1300 


4100 


0.57 ±0.01 


1.70 ± 0.04 


2.41 ±0.07 


0.03 


0.3 


1.5 


0.24 


1300 


2077 


0.57 ±0.01 


1.50 ±0.04 


1.91 ±0.07 



Table 7: (Continued from Table 6) 
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Fig. 10. — Plots for the non-dimensional growth rate Amax (left) and the non-dimensional horizontal 
wavenumber /max (right) of the fastest-growing mode, for Pr= r. Left: We see that Amax scales with 
Pr. The dotted line shows the asymptotic solution of (A7), multiplied by Pr. Right: The horizontal 
wavenumber rapidly becomes independent of Pr. The fastest-g rowing wavclsngth is 27r//max; 

so of 

the order of 13-20d at planetary parameter regimes. The dotted line shows the asymptotic solution 
of (lATl). 
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Fig. 11. — Illustration of the method used to select tstart) applied to the simulation with Pr = r = 
0.1, Rq^ = 1.75. We first find the time t 

min when the total kinetic energy (post saturation), has 
its first local minimum. We then define the width of the peak At as the time interval elapsed since 
the last time the total kinetic energy had the same value. Finally, we define tstart = tmin + 2At. 



- 44 - 



10 



is: 




500 1000 1500 2000 2500 3000 



Fig. 12. — Analysis of the simulation with Pr = r = 0.1, Rq = 1.75. This plot shows the 
temporal evolution of the total kinetic energy in the system (thick black line), half that quantity 
(thin black line), as well as the total kinetic energy in the six highest-amplitude gravity wave mode 
families (see main text for definition and notation). Around t = 2000 the kinetic energy in the 012 
mode family reaches 1/2 the total kinetic energy in the system, at which point it clearly begins to 
dominate the system's transport rates. 



